LIBRARIES

#1.CHM PITFREE 0.5 m

LIDAR_dir <- file.path(system.file("extdata", package = "LadderFuelsR"))
LIDAR_out<-file.path(system.file("extdata", package = "LadderFuelsR"))

las_list <-list.files (LIDAR_dir, pattern = "*.las", full.names = TRUE)
las_list1 <-list.files (LIDAR_dir, pattern = "*.las", full.names = FALSE)
short_name1<-stri_sub(las_list1, 1,-5)

lidar_maps<-lapply(las_list, function(X) lidR::readLAS(X, filter = "-drop_z_below 0"))

#plot(lidar_maps[[1]])

chm_pitfree_list<-list()
outputnames2<-NULL

for (i in seq_along(lidar_maps)) {
  
  lidar_maps1<-lidar_maps[[i]]
  
  chm_pitfree<- grid_canopy(lidar_maps1, res=0.5,pitfree( c(0,2,5,10,15,20,25,30,35,40), c(0,1.5), subcircle=0.15))
  chm_pitfree[chm_pitfree > 40] <- NA
  chm_pitfree[chm_pitfree < 0] <- 0
  chm_pitfree1 <- projectRaster(chm_pitfree, crs=26916)
  
  chm_pitfree_list[[i]]<-chm_pitfree1
  
  outputnames2<- paste0(LIDAR_out ,"/", short_name1[i] ,"_CHM05m.tif",collapse = NULL,sep="")
  writeRaster(chm_pitfree_list[[i]],outputnames2, overwrite=T)
  
}
col <- height.colors(25)
plot(chm_pitfree_list[[1]],col=col)

#2.DETECT TREE TOPS

ttop_dir<- file.path(system.file("extdata", package = "LadderFuelsR"))

las_dir <- file.path(system.file("extdata", package = "LadderFuelsR"))
las_list <-list.files (las_dir, pattern = "*.las", full.names = TRUE)
las_list1 <-list.files (las_dir, pattern = "*.las", full.names = FALSE)
las_files<-lapply(las_list, function(X) lidR::readLAS(X, filter = "-drop_z_below 0"))


names1<-list()
filename_post2<-list()
chm_ttop_list2<-list()

for (j in seq_along(las_list)){
  
  short_name<-stri_sub(las_list1[j], 1, -5)
  names1<-rbind(names1,short_name)
  filename_post2 <- paste0(names1,"_multi2_tree_tops" )
  
  las_file <- lidR::readLAS(las_list[j],"-drop_z_below 0")
  
  # parameters
  ws= 2.5
  hmin = 2
  res=0.5
  ttops_multichm = find_trees(las_file, multichm(res = res, dist_2d = 2,ws= ws, layer_thickness = 0.3,dist_3d = 1, hmin = hmin))
  #ttops_multichm = find_trees(las_file, multichm(res = 2, dist_2d = 2,ws= 5, layer_thickness = 0.3,dist_3d = 1, hmin = 4))
  proj4string(ttops_multichm) <- CRS('+init=EPSG:26916')
  
  chm_ttop_list2[[j]]<-ttops_multichm
  
  writeOGR(chm_ttop_list2[[j]], dsn=ttop_dir ,layer = filename_post2[[j]],driver = 'ESRI Shapefile',
           overwrite_layer = T)
}

las1<-las_files[[1]]
chm1<-chm_pitfree_list[[1]]
ttops1<-chm_ttop_list2[[1]]

#plot(chm1, col = height.colors(50))
#plot(ttops1, add = TRUE, pch= 16, col = "black", main = "Tree tops over the CHM")

# Create an rgl point cloud
x<-add_treetops3d(plot(las_files[[1]], bg = "white", size = 4), ttops1)

# Customize the plot orientation
rgl.viewpoint(theta = 0, phi = 0, fov = 60, zoom = 0.75)

# Convert the rgl scene to an HTML widget
rglwidget(elementId = "x", width = 800, height = 600)

#3.CROWNS (Silva) with locate_trees (ws: MULTICHM FIXED)

crowns_dir <- file.path(system.file("extdata", package = "LadderFuelsR"))

ttop_dir <- file.path(system.file("extdata", package = "LadderFuelsR"))
ttop_list <-list.files (ttop_dir, pattern = "*_multi2_tree_tops.shp", full.names = TRUE)
ttops_multichm<-lapply(ttop_list, function (x) st_read(x))
## Reading layer `Eglin_zone1__crown2m_silva_multi2_tree_tops' from data source 
##   `D:\OLGA\R_library\LadderFuelsR\extdata\Eglin_zone1__crown2m_silva_multi2_tree_tops.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 886 features and 2 fields
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: 525912.8 ymin: 3378523 xmax: 526028.2 ymax: 3378614
## Projected CRS: NAD83 / UTM zone 16N
## Reading layer `Eglin_zone1_000000_multi2_tree_tops' from data source 
##   `D:\OLGA\R_library\LadderFuelsR\extdata\Eglin_zone1_000000_multi2_tree_tops.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 886 features and 2 fields
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: 525912.2 ymin: 3378523 xmax: 526028.2 ymax: 3378614
## Projected CRS: NAD83 / UTM zone 16N
chm_dir <- file.path(system.file("extdata", package = "LadderFuelsR"))
chm_list <-list.files (chm_dir, pattern = "_CHM05m.tif", full.names = TRUE)
las_dir <- file.path(system.file("extdata", package = "LadderFuelsR"))
las_list <-list.files (las_dir, pattern = "*.las", full.names = TRUE)


names1<-list()
filename_post1<-list()
crowns_silva_list<-list()

for (j in seq_along(las_list)){
  
  short_name<-gsub(".*/","",las_list[j])
  short_name1<-stri_sub(short_name, 1, -11)
  names1<-rbind(names1,short_name1)
  filename_post1 <- paste0(names1,"crown2m_silva.las" )

  chm05m <- raster(chm_list[j])
  ttops<-ttops_multichm[[j]]
  las_file <- lidR::readLAS(las_list[j])
  
  algo_silva1 <-silva2016(chm05m, ttops, max_cr_factor = 0.6, exclusion = 0.3, ID = "treeID")
  crowns_silva_las1 <-segment_trees(las_file, algo_silva1, attribute = "treeID", uniqueness = "incremental")
  crowns_silva_las2<-filter_poi(crowns_silva_las1, !is.na(treeID))
  
  crowns_silva_list[[j]]<-crowns_silva_las2

  lidR::writeLAS(crowns_silva_list[[j]], file = paste0(crowns_dir, "/", filename_post1[[j]], sep=""))
}

my_palette <- colorRampPalette(col)
x1<-plot(crowns_silva_list[[1]], color = "treeID", pal = my_palette, bg = "white")

# Customize the plot orientation
rgl.viewpoint(theta = 0, phi = 0, fov = 10, zoom = 0.75)

# Convert the rgl scene to an HTML widget
rglwidget(elementId = "x1", width = 800, height = 600)

FUNCTION CROWN METRICS

custom_crown_metrics <- function(z, i) { # user-defined function
  metrics <- list(
     dz = 1, 
     th = 1,
     z_max = max(z),   # max height
     z_min = min(z),   # min height
     z_mean = mean(z),   # mean height
     z_sd = sd(z), # vertical variability of points
     z_q1=quantile(z, probs = 0.01),
     z_q5=quantile(z, probs = 0.05),
    z_q25=quantile(z, probs = 0.25),
    z_q50=quantile(z, probs = 0.50),
    z_q75=quantile(z, probs = 0.75),
    z_q95=quantile(z, probs = 0.95),
     crr=(mean(z)-min(z))/(max(z)-min(z))
   )
   return(metrics) # output
}
ccm = ~custom_crown_metrics(z = Z, i = Intensity)

#4.TREE METRICS CROWNS leafR (Silva)

stats_out<-file.path(system.file("extdata", package = "LadderFuelsR"))
las_dir <- file.path(system.file("extdata", package = "leafR"))
las_list <-list.files (las_dir, pattern = "*_crown2m_silva.las", full.names = TRUE)

filename_post<-NULL
names2<-NULL
crown_list<-list()

for (j in seq_along(las_list)){
  
  short_name<-stri_sub(las_list[j], 1, -5)
  short_name1<-gsub(".*/","",short_name)
  names2<-rbind(names2,short_name1)
  filename_post <- paste0(names2,"_metrics_leafR.txt" )
  
  normlas_file<-lidR::readLAS(las_list[[j]])
  
  # Setting the xyz coordinates and subsetting the data
  
  xyzi<-subset(normlas_file@data[,c(1:3,5,17)])
  
  xyzi1<-filter(xyzi, Z >= 0.3)
  xyzi1$treeID<-as.factor(xyzi1$treeID)
  TreesMetrics<-CrownMetrics(xyzi1)
  TreesMetrics1<-do.call(rbind, TreesMetrics)
  TreesMetrics3 <- t(TreesMetrics1)
  
  crown_list[[j]]<-TreesMetrics3
  
  write.table(crown_list[[j]], paste0(stats_out, "/", filename_post[[j]] , sep=""), sep="\t", row.names = FALSE)
  
}
## .................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................
metrics1<-crown_list[[1]]
head(metrics1)
##      Tree TotalReturns ETOP     NTOP    EMIN     NMIN     EMAX    NMAX   
## [1,] 1    51           526023.5 3378613 526021.8 526024.9 3378611 3378614
## [2,] 2    136          525997.4 3378587 525995.8 526000   3378585 3378590
## [3,] 3    84           526020.6 3378529 526019.6 526022.4 3378527 3378532
## [4,] 4    100          525976.2 3378550 525974.9 525977.9 3378548 3378552
## [5,] 5    101          525965.9 3378556 525963   525967.6 3378554 3378557
## [6,] 6    325          525998.2 3378545 525997   526000.1 3378543 3378547
##      EWIDTH NWIDTH HMAX  HMIN HMEAN HMEDIAN HMODE HVAR  HSD  HCV   HKUR HSKE 
## [1,] 3.19   2.77   15.77 1.88 11.53 11.31   10.54 8.12  2.85 24.7  4.04 -0.71
## [2,] 4.12   4.66   14.99 3.13 11.28 11.47   14.51 7.85  2.8  24.83 2.35 -0.52
## [3,] 2.85   4.74   14.67 2.86 10.46 11.18   9.53  8.76  2.96 28.31 2.88 -0.88
## [4,] 2.98   4.32   14.66 0.57 10.1  8.84    5.4   13.69 3.7  36.62 2.41 -0.52
## [5,] 4.51   3.42   14.53 0.33 9.49  10.01   0.33  12.58 3.55 37.39 4.01 -1.04
## [6,] 3.09   3.98   14.53 0.33 10.36 11.62   13.62 12.67 3.56 34.34 3.88 -1.12
##      H05TH H10TH H15TH H20TH H25TH H30TH H35TH H40TH H45TH H50TH H55TH H60TH
## [1,] 7.08  8.35  9.3   9.7   9.8   10.17 10.54 10.88 11.25 11.31 11.41 12.52
## [2,] 7.02  7.3   7.71  8.42  9.04  9.92  10.11 10.31 11.22 11.47 12.68 12.84
## [3,] 4.4   5.16  6.66  8.85  9.43  9.53  9.82  10.16 10.64 11.18 11.62 12.12
## [4,] 3.52  5.4   6.37  7.63  7.82  7.91  7.96  8.29  8.48  8.84  12.38 12.87
## [5,] 0.42  6.4   7.11  7.25  7.65  7.87  8.53  8.87  9.56  10.01 10.48 10.7 
## [6,] 0.55  6.83  6.92  7.15  7.28  8.41  9.54  10.19 10.63 11.62 12.2  12.33
##      H65TH H70TH H75TH H80TH H90TH H95TH H99TH IMAX IMIN IMEAN IMEDIAN IMODE
## [1,] 12.68 12.94 13.73 14.91 15.12 15.36 15.77 145  3    51.65 43      3    
## [2,] 13.04 13.35 13.54 14.39 14.51 14.75 14.88 139  0    68.58 74.5    118  
## [3,] 12.27 12.41 12.73 12.98 13.37 14.25 14.54 123  0    42.9  39.5    5    
## [4,] 12.97 13.08 13.56 14    14.27 14.51 14.64 112  0    38.47 33      15   
## [5,] 10.99 11.38 12.33 12.63 13.36 14.25 14.52 137  0    55.96 62      92   
## [6,] 12.77 13.37 13.47 13.64 13.78 13.9  14.35 151  0    55.42 54      7    
##      IVAR    ISD   ICV   IKUR ISKE  I05TH I10TH I15TH I20TH I25TH I30TH I35TH
## [1,] 1464.79 38.27 74.1  1.95 0.35  3.5   7     8     12    17.5  22    28   
## [2,] 1312.94 36.23 52.83 2.09 -0.28 6     12    21.5  33    43.75 49.5  55   
## [3,] 1169.77 34.2  79.72 2.05 0.46  2     3.3   5     6.2   9.5   15.8  21.1 
## [4,] 782.6   27.97 72.72 2.32 0.54  3.95  5.9   9.7   12.8  14.75 17    19   
## [5,] 1387.94 37.26 66.57 1.8  0.06  2     6     11    14    21    29    35   
## [6,] 1420.3  37.69 68.01 1.93 0.25  4     7     11    15    21    27    32   
##      I40TH I45TH I50TH I55TH I60TH I65TH I70TH I75TH I80TH  I90TH I95TH  I99TH 
## [1,] 30    33    43    54    69    71.5  75    86    93.5   100   109.5  129   
## [2,] 62    68.75 74.5  78    84    89    92    95.25 105.75 116   121.25 132.2 
## [3,] 27.6  34.4  39.5  43    48.6  51.95 55    69.75 88     92.7  98.55  117.19
## [4,] 21    29.1  33    39    44.4  49    54.3  62.25 69.15  79.1  84.2   109.03
## [5,] 39    51    62    68    73    76    81    84    97     101   109    131   
## [6,] 38    44.8  54    60    68    74    80    85    101    107.6 117.8  132.52

#5.TREE METRICS CROWNS standard and customized (Silva)

stats_out<-file.path(system.file("extdata", package = "LadderFuelsR"))
las_dir <- file.path(system.file("extdata", package = "LadderFuelsR"))
las_list <-list.files (las_dir, pattern = "*_crown2m_silva.las", full.names = TRUE)

names2<-list()
filename_post4<-list()
crowns_metrics_list<-list()

for (j in seq_along(las_list)){
  
  short_name<-stri_sub(las_list[j], 1, -5)
  short_name1<-gsub(".*/","",short_name)
  names2<-rbind(names2,short_name1)
  filename_post4 <- paste0(names2,"_metrics_convex.shp" )
  
  las_file1 <- lidR::readLAS(las_list[j])
  las_file2<-filter_poi(las_file1, !is.na(treeID))
  las_file3<-filter_poi(las_file2, Z >= 1)
  
  metrics1 = crown_metrics(las_file3,func = .stdtreemetrics, geom = "convex")
  crown_diam<-data.frame(sqrt(metrics1$convhull_area/ pi) * 2)
  names(crown_diam)<-"crown_diam"
  metrics2 = crown_metrics(las_file3,func = ccm, geom = "convex") #concave
  metrics_all <- dplyr::bind_cols(list(metrics1,crown_diam,metrics2))
  metrics_all1 <- metrics_all[,c(1:4,6,10:21)]
  names(metrics_all1)<-c("treeID", "Z", "npoints", "convhull_area", "crown_diam", "z_max", "z_min", "z_mean","z_sd", "z_q1","z_q5", "z_q25","z_q50","z_q75", "z_q95", "crr", "geometry" )
  
  metrics_all2<- st_as_sf(metrics_all1)
  
  crowns_metrics_list[[j]]<-metrics_all2
  
  st_write(crowns_metrics_list[[j]], paste0(dsn=stats_out, "/", filename_post4[[j]],sep=""), driver = 'ESRI Shapefile',append=F)
  
}
## Writing layer `Eglin_zone1__crown2m_crown2m_silva_metrics_convex' to data source 
##   `D:/OLGA/R_library/LadderFuelsR/extdata/Eglin_zone1__crown2m_crown2m_silva_metrics_convex.shp' using driver `ESRI Shapefile'
## Writing 878 features with 16 fields and geometry type Polygon.
## Deleting layer `Eglin_zone1__crown2m_silva_metrics_convex' using driver `ESRI Shapefile'
## Writing layer `Eglin_zone1__crown2m_silva_metrics_convex' to data source 
##   `D:/OLGA/R_library/LadderFuelsR/extdata/Eglin_zone1__crown2m_silva_metrics_convex.shp' using driver `ESRI Shapefile'
## Writing 877 features with 16 fields and geometry type Polygon.
ttops1<-st_as_sf(chm_ttop_list2[[1]])
crowns1<-st_as_sf(crowns_metrics_list[[1]])
ttops_within_crowns <- st_intersection(ttops1, crowns1)

# Set the size of the plotting device
par(mfrow = c(1, 1), mar = c(1, 1, 1, 1), pin = c(5, 4)) 
plot(st_geometry(crowns1), pch = 16, col = "green")
plot(ttops_within_crowns, add = TRUE, pch= 16, col = "darkblue", main = "Tree tops over the crowns")

#6.NO OVERLAPPING CROWN POLYGONS

dir_out<-file.path(system.file("extdata", package = "LadderFuelsR"))
tree_dir <- file.path(system.file("extdata", package = "LadderFuelsR"))
tree_list <-list.files (tree_dir, pattern = "*_crown2m_silva_metrics_convex.shp", full.names = TRUE)
tree_files <-lapply(tree_list, function (X) st_read(X))
## Reading layer `Eglin_zone1__crown2m_crown2m_silva_metrics_convex' from data source `D:\OLGA\R_library\LadderFuelsR\extdata\Eglin_zone1__crown2m_crown2m_silva_metrics_convex.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 878 features and 16 fields
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: 525912.5 ymin: 3378524 xmax: 526028 ymax: 3378614
## Projected CRS: NAD83 / UTM zone 16N
## Reading layer `Eglin_zone1__crown2m_silva_metrics_convex' from data source 
##   `D:\OLGA\R_library\LadderFuelsR\extdata\Eglin_zone1__crown2m_silva_metrics_convex.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 877 features and 16 fields
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: 525912.5 ymin: 3378523 xmax: 526028 ymax: 3378614
## Projected CRS: NAD83 / UTM zone 16N
polys1<-tree_files[[1]]

for(i in 1:nrow(polys1)) {
  other_polys <- st_union(polys1[-i, ]$geometry)
  polys1[i, ]$geometry <- st_difference(polys1[i, ]$geometry, other_polys)
}
st_write(polys1,paste0(tree_dir, "/", "crown2m_silva_convex_no_overlap.shp"), append=FALSE)
## Deleting layer `crown2m_silva_convex_no_overlap' using driver `ESRI Shapefile'
## Writing layer `crown2m_silva_convex_no_overlap' to data source 
##   `D:/OLGA/R_library/LadderFuelsR/extdata/crown2m_silva_convex_no_overlap.shp' using driver `ESRI Shapefile'
## Writing 878 features with 16 fields and geometry type Polygon.
# Set the size of the plotting device
par(mfrow = c(1, 1), mar = c(1, 1, 1, 1), pin = c(5, 4)) 
plot(st_geometry(polys1), pch = 16, col = "orange")
plot(ttops_within_crowns, add = TRUE, pch= 16, col = "blue", main = "Tree tops over the crowns")

#7.CROP LAS FILES WITH CROWN POLYGONS

dir_out<-file.path(system.file("extdata", package = "LadderFuelsR"))

las_dir <- file.path(system.file("extdata", package = "LadderFuelsR"))
las_list <-list.files (las_dir, pattern = "*_crown2m_silva.las", full.names = TRUE)
las_files <-lapply(las_list, function (X) lidR::readLAS(X))

tree_dir <- file.path(system.file("extdata", package = "LadderFuelsR"))
tree_list <-list.files (tree_dir, pattern = "*_no_overlap.shp", full.names = TRUE)
tree_files <-lapply(tree_list, function (X) st_read(X))
## Reading layer `crown2m_silva_convex_no_overlap' from data source 
##   `D:\OLGA\R_library\LadderFuelsR\extdata\crown2m_silva_convex_no_overlap.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 878 features and 16 fields
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: 525912.5 ymin: 3378524 xmax: 526028 ymax: 3378614
## Projected CRS: NAD83 / UTM zone 16N
###################################

names1 <- list()
filename_post3 <- list()

for (j in seq_along(las_list)) {
  short_name <- stri_sub(las_list[j])
  short_name1 <- gsub(".*/", "", short_name)
  short_name2 <- stri_sub(short_name1, 1, 11)
  
  print(short_name2)
  
  if (is.na(short_name2)) {
    # Handle NA values or skip them
    next
  }
  
  names1[[j]] <- short_name2
  filename_post3[[j]] <- paste0(short_name2, "_CROWN.las")
  
  las1 <- las_files[[j]] 
  trees1 <- tree_files[[j]] 
  
  trees_ID <- trees1 %>% dplyr::select(treeID)
  n <- nrow(trees_ID)
  
  crown_cort <- vector("list", length=n)
  
  for (i in 1:n) {
    kk <- trees_ID[i,]
    crown_cort[[i]] = clip_roi(las1, kk)
  
    lidR::writeLAS(crown_cort[[i]], paste0(dir_out,"/", i, "_", filename_post3[j]))
  }
}
## [1] "Eglin_zone1"
my_palette <- colorRampPalette(col)
x2<-plot(crown_cort[[1]], color = "Z", pal = my_palette, bg = "black", size = 2.5)

# Customize the plot orientation
rgl.viewpoint(theta = 0, phi = 0, fov = 60, zoom = 0.75)

# Convert the rgl scene to an HTML widget
rglwidget(elementId = "x2", width = 400, height = 600)

#8.LAI-LAD METRICS BY TREE

stats_out<-file.path(system.file("extdata", package = "LadderFuelsR"))

las_dir1 <- file.path(system.file("extdata", package = "leafR"))
las_list1 <- list.files(las_dir1, pattern = "*_CROWN.las", full.names = TRUE, ignore.case = TRUE)

# create a vector to hold the file names of .las files with more than 3 points
files_with_more_than_3_points <- c()

# loop through each file
for (file in las_list1) {
  las_data <- lidR::readLAS(file)
  las_data1<-filter_poi(las_data, Z >= 1)
  
  # skip to next file if there was a problem reading
  if (is.null(las_data1)) next
  
  # check if it contains more than three points
  if (las_data1@header$`Number of point records` > 3) {
    files_with_more_than_3_points <- c(files_with_more_than_3_points, file)
  }
}

# Creates a data frame of the 3D voxels information (xyz) with Leaf Area Density values
short_name1<-NULL
profile_list<-NULL
lidar_lai_list<-NULL
understory_lai_list<-NULL
LAHV_metric_list<-NULL

for (j in seq_along(files_with_more_than_3_points)){
  
  short_name<-stri_sub(files_with_more_than_3_points[j], 1, -5)
  short_name1<-gsub(".*/","",short_name)
  
  normlas_file<-files_with_more_than_3_points[[j]]

  VOXELS_LAD = lad.voxels(normlas_file, grain.size = 2)
  
  lad_profile = lad.profile(VOXELS_LAD, relative = F)
  lai_tot = lai(lad_profile)
  understory_lai <- lai(lad_profile, min = 0.3, max = 2.5)
  LAHV_metric<- LAHV(lad_profile, LAI.weighting = FALSE, height.weighting = FALSE) # Leaf Area height volume (wood volume, AGB)
  
  lad_profile1 = data.frame(lad_profile, treeID = short_name1)
  lai_tot1 = data.frame(lai_tot, treeID = short_name1)
  understory_lai1 = data.frame(understory_lai, treeID = short_name1)
  LAHV_metric1 = data.frame(LAHV_metric, treeID = short_name1)
  
  profile_list<-rbind(profile_list, lad_profile1)
  lidar_lai_list<-rbind(lidar_lai_list,lai_tot1)
  understory_lai_list <-rbind(understory_lai_list,understory_lai1)
  LAHV_metric_list<-rbind(LAHV_metric_list,LAHV_metric1)
}

write.table(profile_list,file =paste(stats_out, "/","alltrees_LAD_profile_voxels2m",".txt", sep=""), sep="\t", row.names = FALSE)
write.table(lidar_lai_list,file =paste(stats_out, "/","alltrees_LAI_profile_voxels2m_",".txt", sep=""), sep="\t", row.names = FALSE)
write.table(understory_lai_list,file =paste(stats_out, "/","alltrees_understory_profile_voxels2m",".txt", sep=""), sep="\t", row.names = FALSE)
write.table(LAHV_metric_list,file =paste(stats_out, "/","alltrees_LAHV_profile_voxels2m",".txt", sep=""), sep="\t", row.names = FALSE)
head(profile_list)
##   height        lad              treeID
## 1    1.5 0.01724822 1_Eglin_zone1_CROWN
## 2    2.5 0.00000000 1_Eglin_zone1_CROWN
## 3    3.5 0.00000000 1_Eglin_zone1_CROWN
## 4    4.5 0.00000000 1_Eglin_zone1_CROWN
## 5    5.5 0.00000000 1_Eglin_zone1_CROWN
## 6    6.5 0.10136628 1_Eglin_zone1_CROWN

#9.DEPURATING TREE LAD PROFILES (> 3 HEIGHT VALUES)

LAD_path<- file.path(system.file("extdata", package = "LadderFuelsR"),"alltrees_LAD_profile_voxels2m.txt")
stats_tot<- read.table(LAD_path,sep="\t", header=T)
cols <- c('treeID')
stats_tot[cols] <- lapply(stats_tot[cols], function (x) as.factor(x))
stats_tot$lad<-round(stats_tot$lad,digits = 4)

cases <- data.frame(table(stats_tot$treeID))
cases1 <-cases[cases$Freq > 3, ]
names(cases1)<-c("treeID", "Freq")

stats_tot1 <- stats_tot[stats_tot$treeID %in% cases1$treeID, ]
stats_tot2 <- data.frame(stats_tot1 %>% replace(is.na(.), 0.01))

output_file <- file.path(system.file("extdata", package = "LadderFuelsR"), "LAD_profiles.txt")
write.table(stats_tot2, file = output_file, sep = "\t", row.names = FALSE)

#10.GAPS AND FUEL LAYERS BASE HEIGHT (FBH)

# LAD profiles derived from normalized ALS data after applying [lad.profile()] function
data_path <- file.path(system.file("extdata", package = "LadderFuelsR"), "LAD_profiles.txt")
LAD_profiles <- read.table(data_path, sep = "\t", header = TRUE)
LAD_profiles$treeID <- factor(LAD_profiles$treeID)

trees_name1 <- as.character(LAD_profiles$treeID)
trees_name2 <- factor(unique(trees_name1))

metrics_precentile_list<-list()
for (i in levels(trees_name2)) {
  tree2 <- LAD_profiles |> dplyr::filter(treeID == i)
  metrics_precentil <- get_gaps_fbhs(tree2)
  metrics_precentile_list[[i]] <- metrics_precentil
}

metrics_all_percentil <- dplyr::bind_rows(metrics_precentile_list)
metrics_all_percentil$treeID <- factor(metrics_all_percentil$treeID)

# Remove the row with all NA values from the original data frame
# First remove "treeID" and "treeID1" columns
tree_metrics_no_treeID <- metrics_all_percentil[, -which(names(metrics_all_percentil) == c("treeID","treeID1"))]

# Check if any row has all NA values
rows_with_all_NA_or_zero <- apply(tree_metrics_no_treeID, 1, function(row) all(is.na(row) | row == 0))

# Get the row index with all NA values
row_index <- which(rows_with_all_NA_or_zero)

# Remove the row with all NA values from the original data frame
if (length(row_index) > 0) {
  tree_metrics_filtered <- metrics_all_percentil[-row_index, ]
} else {
  tree_metrics_filtered <- metrics_all_percentil
}
write.table(tree_metrics_filtered, file= file.path(system.file("extdata", package = "LadderFuelsR"), "1_gaps_fbhs_metrics.txt"),
sep = "\t",row.names = FALSE)

head(tree_metrics_filtered)
##                           treeID cbh1 gap1 gap2 cbh2 cbh3 cbh4 gap_lad1
## height...1   1_Eglin_zone1_CROWN  1.5  2.5  5.5  6.5  9.5 14.5        0
## height...2  10_Eglin_zone1_CROWN  5.5  1.5  4.5 10.5 13.5 <NA>        0
## height...3 100_Eglin_zone1_CROWN  1.5  2.5  4.5  5.5  8.5 11.5        0
## height...4 101_Eglin_zone1_CROWN  4.5  1.5  3.5  6.5 11.5 <NA>        0
## height...5 102_Eglin_zone1_CROWN  5.5  1.5  4.5  8.5 11.5 <NA>        0
## height...6 103_Eglin_zone1_CROWN  1.5  2.5  4.5  5.5  7.5 10.5        0
##            gap_lad2 cbh_perc1 cbh_perc2 cbh_perc3 cbh_perc4 cbh_lad1 cbh_lad2
## height...1        0        50        75        90        50   0.0172   0.1014
## height...2        0        75        75       100        NA   0.2310   0.0608
## height...3        0        50        50        95        75   0.0239   0.0143
## height...4        0        75       100        50        NA   0.1032   0.2544
## height...5        0        90       100        50        NA   0.1754   0.2338
## height...6        0        50        75        95       100   0.0230   0.0445
##            cbh_lad3 cbh_lad4 max_height treeID1 gap3 gap_lad3 cbh5 cbh_perc5
## height...1   0.1933   0.0333       15.5       1 <NA>       NA <NA>        NA
## height...2   0.6150       NA       14.5      10  9.5        0 <NA>        NA
## height...3   0.1705   0.0924       12.5     100 <NA>       NA <NA>        NA
## height...4   0.0668       NA       12.5     101 <NA>       NA <NA>        NA
## height...5   0.0274       NA       12.5     102 <NA>       NA <NA>        NA
## height...6   0.1744   0.3261       12.5     103 <NA>       NA 11.5        50
##            cbh_lad5 cbh6 cbh_perc6 cbh_lad6 gap4 gap5 gap_lad4 gap_lad5
## height...1       NA <NA>        NA       NA <NA> <NA>       NA       NA
## height...2       NA <NA>        NA       NA <NA> <NA>       NA       NA
## height...3       NA <NA>        NA       NA <NA> <NA>       NA       NA
## height...4       NA <NA>        NA       NA <NA> <NA>       NA       NA
## height...5       NA <NA>        NA       NA <NA> <NA>       NA       NA
## height...6   0.0343 <NA>        NA       NA <NA> <NA>       NA       NA

#11.DISTANCE BETWEEN FUEL LAYERS

# Tree metrics derived from get_gaps_fbhs() function
gaps_fbhs_metrics_path <- system.file("extdata", "1_gaps_fbhs_metrics.txt", package = "LadderFuelsR")
gaps_fbhs_metrics <- read.table(gaps_fbhs_metrics_path, sep="\t", header=TRUE)

gaps_fbhs_metrics$treeID <- factor(gaps_fbhs_metrics$treeID)
trees_name1 <- as.character(gaps_fbhs_metrics$treeID)
trees_name2 <- factor(unique(trees_name1))

metrics_distance_list <- list()

for (i in levels(trees_name2)) {

  # Filter data for each tree
  tree2 <- gaps_fbhs_metrics |> dplyr::filter(treeID == i)

  # Get distance metrics for each tree
  metrics_distance <- get_distance(tree2)
  metrics_distance_list[[i]] <- metrics_distance
}

# Combine the individual data frames
metrics_all_distance <- dplyr::bind_rows(metrics_distance_list)

distance_path <- file.path(system.file("extdata", package = "LadderFuelsR"), "2_distance_metrics.txt")
write.table(metrics_all_distance, file=distance_path, sep="\t", row.names=FALSE)

metrics_all_distance <- metrics_all_distance[, order(names(metrics_all_distance))]
head(metrics_all_distance)
##   cbh1 cbh2 cbh3 cbh4 cbh5 cbh6 dist1 dist2 dist3 gap1 gap2 gap3 gap4 gap5
## 1  1.5  6.5  9.5 14.5   NA   NA     4    NA    NA  2.5  5.5   NA   NA   NA
## 2  5.5 10.5 13.5   NA   NA   NA     4     1    NA  1.5  4.5  9.5   NA   NA
## 3  1.5  5.5  8.5 11.5   NA   NA     3    NA    NA  2.5  4.5   NA   NA   NA
## 4  4.5  6.5 11.5   NA   NA   NA     3    NA    NA  1.5  3.5   NA   NA   NA
## 5  5.5  8.5 11.5   NA   NA   NA     4    NA    NA  1.5  4.5   NA   NA   NA
## 6  1.5  5.5  7.5 10.5 11.5   NA     3    NA    NA  2.5  4.5   NA   NA   NA
##   Hdist1 Hdist2 Hdist3 max_height                treeID treeID1
## 1    5.5     NA     NA       15.5   1_Eglin_zone1_CROWN       1
## 2    4.5    9.5     NA       14.5  10_Eglin_zone1_CROWN      10
## 3    4.5     NA     NA       12.5 100_Eglin_zone1_CROWN     100
## 4    3.5     NA     NA       12.5 101_Eglin_zone1_CROWN     101
## 5    4.5     NA     NA       12.5 102_Eglin_zone1_CROWN     102
## 6    4.5     NA     NA       12.5 103_Eglin_zone1_CROWN     103

#12.FUEL LAYERS DEPTH

library(dplyr)
library(magrittr)

# LAD profiles derived from normalized ALS data after applying [lad.profile()] function
data_path <- file.path(system.file("extdata", package = "LadderFuelsR"), "LAD_profiles.txt")
LAD_profiles <- read.table(data_path, sep = "\t", header = TRUE)
LAD_profiles$treeID <- factor(LAD_profiles$treeID)

# Tree metrics derived from get_distance() function
distance_path <- file.path(system.file("extdata", package = "LadderFuelsR"), "2_distance_metrics.txt")
distance_metrics <- read.table(distance_path, sep = "\t", header = TRUE)
distance_metrics$treeID <- factor(distance_metrics$treeID)

metrics_depth_list <- list()

for (i in levels(LAD_profiles$treeID)){

  tree1 <- LAD_profiles |> dplyr::filter(treeID == i)
  tree2 <- distance_metrics |> dplyr::filter(treeID == i)

  # Get depths for each tree
  metrics_depth <- get_depths(tree1, tree2)
  metrics_depth_list[[i]] <- metrics_depth
}

# Combine the individual data frames
metrics_all_depth <- dplyr::bind_rows(metrics_depth_list)

depth_path <- file.path(system.file("extdata", package = "LadderFuelsR"), "3_depth_metrics.txt")
write.table(metrics_all_depth, file = depth_path, sep = "\t", row.names = FALSE)

metrics_all_depth <- metrics_all_depth[, order(names(metrics_all_depth))]
head(metrics_all_depth)
##   cbh1 cbh2 cbh3 cbh4 cbh5 cbh6 depth0 depth1 depth2 depth3 dist1 dist2 dist3
## 1  1.5  6.5  9.5 14.5   NA   NA      1      8     NA     NA     4    NA    NA
## 2  5.5 10.5 13.5   NA   NA   NA      0      4      3     NA     4     1    NA
## 3  1.5  5.5  8.5 11.5   NA   NA      1      6     NA     NA     3    NA    NA
## 4  4.5  6.5 11.5   NA   NA   NA      0      7     NA     NA     3    NA    NA
## 5  5.5  8.5 11.5   NA   NA   NA      0      6     NA     NA     4    NA    NA
## 6  1.5  5.5  7.5 10.5 11.5   NA      1      6     NA     NA     3    NA    NA
##   gap1 gap2 gap3 gap4 gap5 Hdepth0 Hdepth1 Hdepth2 Hdepth3 Hdist1 Hdist2 Hdist3
## 1  2.5  5.5   NA   NA   NA     1.5    14.5      NA      NA    5.5     NA     NA
## 2  1.5  4.5  9.5   NA   NA     0.5     8.5    13.5      NA    4.5    9.5     NA
## 3  2.5  4.5   NA   NA   NA     1.5    11.5      NA      NA    4.5     NA     NA
## 4  1.5  3.5   NA   NA   NA     0.5    11.5      NA      NA    3.5     NA     NA
## 5  1.5  4.5   NA   NA   NA     0.5    11.5      NA      NA    4.5     NA     NA
## 6  2.5  4.5   NA   NA   NA     1.5    11.5      NA      NA    4.5     NA     NA
##   max_height                treeID treeID1
## 1       15.5   1_Eglin_zone1_CROWN       1
## 2       14.5  10_Eglin_zone1_CROWN      10
## 3       12.5 100_Eglin_zone1_CROWN     100
## 4       12.5 101_Eglin_zone1_CROWN     101
## 5       12.5 102_Eglin_zone1_CROWN     102
## 6       12.5 103_Eglin_zone1_CROWN     103

#13.PLOTS GAPS AND FUEL LAYERS BASE HEIGHT (FBH)

library(LadderFuelsR)
library(ggplot2)
library(lattice)

# LAD profiles derived from normalized ALS data after applying [lad.profile()] function
data_path <- file.path(system.file("extdata", package = "LadderFuelsR"), "LAD_profiles.txt")
LAD_profiles <- read.table(data_path, sep = "\t", header = TRUE)
LAD_profiles$treeID <- factor(LAD_profiles$treeID)

# Tree metrics derived from get_depths() function
depth_path <- file.path(system.file("extdata", package = "LadderFuelsR"), "3_depth_metrics.txt")
depth_metrics <- read.table(depth_path, sep = "\t", header = TRUE)
depth_metrics$treeID <- factor(depth_metrics$treeID)

# Generate plots for gaps and fbhs
plots_gaps_fbhs <- get_plots_gap_fbh(LAD_profiles, depth_metrics)

# Save plots for each tree
for (name in names(plots_gaps_fbhs)) {
  #print(plots_gaps_fbhs[[name]])  # print the plot
  ggsave(file.path(system.file("extdata", package = "LadderFuelsR"),  paste0( name, "_gap_fbh", ".tiff")), plot = plots_gaps_fbhs[[name]])
}

par(mfrow = c(3, 1))
# Plot the first panel
plot(plots_gaps_fbhs[[1]],width = 5, height = 5)

# Plot the second panel
plot(plots_gaps_fbhs[[2]],width = 5, height = 5)

# Plot the third panel
plot(plots_gaps_fbhs[[3]],width = 5, height = 5)

#14.FUEL LAYERS BASE HEIGHT (FBH) AFTER REMOVING DISTANCES = 1

library(SSBtools)
library(dplyr)
library(magrittr)

# Tree metrics derived from get_depths() function
depth_path <- file.path(system.file("extdata", package = "LadderFuelsR"), "3_depth_metrics.txt")
depth_metrics <- read.table(depth_path, sep = "\t", header = TRUE)
depth_metrics$treeID <- factor(depth_metrics$treeID)

trees_name1 <- as.character(depth_metrics$treeID)
trees_name2 <- factor(unique(trees_name1))

fbh_corr_list <- list()

for (i in levels(trees_name2)){

  # Filter data for each tree
  tree3 <- depth_metrics |> dplyr::filter(treeID == i)

  # Get real fbh for each tree
  fbh_corr <- get_real_fbh(tree3)

  # Store fbh values in a list
  fbh_corr_list[[i]] <- fbh_corr
}

# Combine fbh values for all trees
fbh_corr_all <- dplyr::bind_rows(fbh_corr_list)
fbh_corr_all$treeID <- factor(fbh_corr_all$treeID)

# Reorder columns
# Get original column names
original_column_names <- colnames(fbh_corr_all)

# Specify prefixes
prefixes <- c("treeID", "Hdist", "Hcbh", "Hdepth", "dist", "depth", "max_height")

# Initialize vector to store new order
new_order <- c()

# Loop over prefixes
for (prefix in prefixes) {
  # Find column names matching the current prefix
  matching_columns <- grep(paste0("^", prefix), original_column_names, value = TRUE)
  # Append to the new order
  new_order <- c(new_order, matching_columns)
}

# Reorder values
fbh_corr_all <- fbh_corr_all[, new_order]

# Save the reordered data
fbh_corr_path <- file.path(system.file("extdata", package = "LadderFuelsR"), "4_fbh_metrics_corr.txt")
write.table(fbh_corr_all, file = fbh_corr_path, sep = "\t", row.names = FALSE)
head(fbh_corr_all)
##                             treeID treeID1 Hdist1 Hdist2 Hdist3 Hcbh1 Hcbh2
## new_Hcbh...1   1_Eglin_zone1_CROWN       1    5.5     NA     NA   1.5   6.5
## new_Hcbh...2  10_Eglin_zone1_CROWN      10    4.5    9.5     NA   5.5   5.5
## new_Hcbh...3 100_Eglin_zone1_CROWN     100    4.5     NA     NA   1.5   5.5
## new_Hcbh...4 101_Eglin_zone1_CROWN     101    3.5     NA     NA   4.5   4.5
## new_Hcbh...5 102_Eglin_zone1_CROWN     102    4.5     NA     NA   5.5   5.5
## new_Hcbh...6 103_Eglin_zone1_CROWN     103    4.5     NA     NA   1.5   5.5
##              Hcbh3 Hcbh4 Hcbh5 Hcbh6 Hdepth1 Hdepth2 Hdepth3 Hdepth4 dist1
## new_Hcbh...1   6.5   6.5    NA    NA     1.5    14.5      NA      NA     4
## new_Hcbh...2   5.5    NA    NA    NA     8.5    13.5      NA      NA     4
## new_Hcbh...3   5.5   5.5    NA    NA     1.5    11.5      NA      NA     3
## new_Hcbh...4   4.5    NA    NA    NA    11.5      NA      NA      NA     3
## new_Hcbh...5   5.5    NA    NA    NA    11.5      NA      NA      NA     4
## new_Hcbh...6   5.5   5.5   5.5    NA     1.5    11.5      NA      NA     3
##              dist2 dist3 depth1 depth2 depth3 depth4 max_height
## new_Hcbh...1    NA    NA      1      8     NA     NA       15.5
## new_Hcbh...2     1    NA      4      3     NA     NA       14.5
## new_Hcbh...3    NA    NA      1      6     NA     NA       12.5
## new_Hcbh...4    NA    NA      7     NA     NA     NA       12.5
## new_Hcbh...5    NA    NA      6     NA     NA     NA       12.5
## new_Hcbh...6    NA    NA      1      6     NA     NA       12.5

#15.FUEL LAYERS DEPTH AFTER REMOVING DISTANCES = 1

library(dplyr)
library(magrittr)
library(tidyr)
# Tree metrics derived from get_real_fbh() function
fbhcor_path <- file.path(system.file("extdata", package = "LadderFuelsR"), "4_fbh_metrics_corr.txt")
fbh_metrics_corr <- read.table(fbhcor_path, sep = "\t", header = TRUE)
fbh_metrics_corr$treeID <- factor(fbh_metrics_corr$treeID)

trees_name1 <- as.character(fbh_metrics_corr$treeID)
trees_name2 <- factor(unique(trees_name1))

depth_metrics_corr <- lapply(levels(trees_name2), function(i) {
  # Filter data for each tree
  tree2 <- fbh_metrics_corr |> dplyr::filter(treeID == i)
  # Get real depths for each tree
  get_real_depths(tree2)
})

# Combine depth values for all trees
depth_corr_all <- dplyr::bind_rows(depth_metrics_corr)
depth_corr_path <- file.path(system.file("extdata", package = "LadderFuelsR"), "5_tree_depth_metrics_corr.txt")
write.table(depth_corr_all, file = depth_corr_path, sep = "\t", row.names = FALSE)

head(depth_corr_all)
##   treeID1                treeID Hcbh1 Hcbh2 dptf1 dptf2 Hdptf1 Hdptf2 dist1
## 1       1   1_Eglin_zone1_CROWN   1.5   6.5     1     8    1.5   14.5     4
## 2      10  10_Eglin_zone1_CROWN   5.5    NA     8    NA   13.5     NA     4
## 3     100 100_Eglin_zone1_CROWN   1.5   5.5     1     6    1.5   11.5     3
## 4     101 101_Eglin_zone1_CROWN   4.5    NA     7    NA   11.5     NA     3
## 5     102 102_Eglin_zone1_CROWN   5.5    NA     6    NA   11.5     NA     4
## 6     103 103_Eglin_zone1_CROWN   1.5   5.5     1     6    1.5   11.5     3
##   Hdist1 max_height Hcbh3 dptf3 Hdptf3 dist2 dist3 Hdist2 Hdist3
## 1    5.5       15.5    NA    NA     NA    NA    NA     NA     NA
## 2    4.5       14.5    NA    NA     NA    NA    NA     NA     NA
## 3    4.5       12.5    NA    NA     NA    NA    NA     NA     NA
## 4    3.5       12.5    NA    NA     NA    NA    NA     NA     NA
## 5    4.5       12.5    NA    NA     NA    NA    NA     NA     NA
## 6    4.5       12.5    NA    NA     NA    NA    NA     NA     NA

#16.FUEL LAYERS DISTANCES (> 1 M) AND CBH BASE ON MAXIMUM- AND LAST- DISTANCE

library(dplyr)
library(magrittr)
library(stringr)

# Tree metrics derived from get_real_depths() function
depth_corr_path <- file.path(system.file("extdata", package = "LadderFuelsR"), "5_tree_depth_metrics_corr.txt")
depth_metrics_corr <- read.table(depth_corr_path, sep = "\t", header = TRUE)
depth_metrics_corr$treeID <- factor(depth_metrics_corr$treeID)

trees_name1 <- as.character(depth_metrics_corr$treeID)
trees_name2 <- factor(unique(trees_name1))

distance_metrics_corr <- lapply(levels(trees_name2), function(i) {
  # Filter data for each tree
  tree2 <- depth_metrics_corr |> dplyr::filter(treeID == i)
  # Get effective gap for each tree
  get_effective_gap(tree2)
})

# Combine the individual data frames
distances_corr_all <- dplyr::bind_rows(distance_metrics_corr)

# =======================================================================#
# REORDER COLUMNS:
# =======================================================================#
# Get original column names
original_column_names <- colnames(distances_corr_all)

# Specify prefixes
prefixes <- c("treeID", "Hcbh", "dptf", "Hdptf", "effdist", "dist", "Hdist", "max_Hcbh", "max_dptf", "max_Hdptf", "last_Hcbh", "last_dptf", "last_Hdptf", "max_height")

# Initialize vector to store new order
new_order <- c()

# Loop over prefixes
for (prefix in prefixes) {
  # Find column names matching the current prefix
  matching_columns <- grep(paste0("^", prefix), original_column_names, value = TRUE)

  # Extract numeric suffixes and order the columns based on these suffixes
  numeric_suffixes <- as.numeric(gsub(paste0("^", prefix), "", matching_columns))
  matching_columns <- matching_columns[order(numeric_suffixes)]

  # Append to new order
  new_order <- c(new_order, matching_columns)
}

# Reorder values
distances_corr_all1 <- distances_corr_all[, new_order]
distance_corr_path <- file.path(system.file("extdata", package = "LadderFuelsR"), "6_tree_distances_metrics_corr.txt")
write.table(distances_corr_all1, file = distance_corr_path, sep = "\t", row.names = FALSE)

head(distances_corr_all1)
##   treeID1                treeID Hcbh1 Hcbh2 Hcbh3 dptf1 dptf2 dptf3 Hdptf1
## 1       1   1_Eglin_zone1_CROWN   1.5   6.5    NA     1     8    NA    1.5
## 2      10  10_Eglin_zone1_CROWN   5.5    NA    NA     8    NA    NA   13.5
## 3     100 100_Eglin_zone1_CROWN   1.5   5.5    NA     1     6    NA    1.5
## 4     101 101_Eglin_zone1_CROWN   4.5    NA    NA     7    NA    NA   11.5
## 5     102 102_Eglin_zone1_CROWN   5.5    NA    NA     6    NA    NA   11.5
## 6     103 103_Eglin_zone1_CROWN   1.5   5.5    NA     1     6    NA    1.5
##   Hdptf2 Hdptf3 effdist1 effdist2 dist1 dist2 dist3 Hdist1 Hdist2 Hdist3
## 1   14.5     NA        4       NA     4    NA    NA    5.5     NA     NA
## 2     NA     NA        4       NA     4    NA    NA    4.5     NA     NA
## 3   11.5     NA        3       NA     3    NA    NA    4.5     NA     NA
## 4     NA     NA        3       NA     3    NA    NA    3.5     NA     NA
## 5     NA     NA        4       NA     4    NA    NA    4.5     NA     NA
## 6   11.5     NA        3       NA     3    NA    NA    4.5     NA     NA
##   max_Hcbh max_dptf max_Hdptf last_Hcbh last_dptf last_Hdptf max_height
## 1      6.5        8      14.5       6.5         8       14.5       15.5
## 2      5.5        8      13.5       5.5         8       13.5       14.5
## 3      5.5        6      11.5       5.5         6       11.5       12.5
## 4      4.5        7      11.5       4.5         7       11.5       12.5
## 5      5.5        6      11.5       5.5         6       11.5       12.5
## 6      5.5        6      11.5       5.5         6       11.5       12.5

#17.FUELS LAD PERCENTAGE (> 25 %) AND CBH BASED ON MAXIMUM LAD PERCENTAGE

library(dplyr)
library(magrittr)

# LAD profiles derived from normalized ALS data after applying [lad.profile()] function
data_path <- system.file("extdata", "LAD_profiles.txt", package = "LadderFuelsR")
LAD_profiles <- read.table(data_path, sep = "\t", header = TRUE)
LAD_profiles$treeID <- factor(LAD_profiles$treeID)

# Tree metrics derived from get_effective_gap() function
distcor_path <- file.path(system.file("extdata", package = "LadderFuelsR"), "6_tree_distances_metrics_corr.txt")
distances_metrics_corr <- read.table(distcor_path, sep = "\t", header = TRUE)
distances_metrics_corr$treeID <- factor(distances_metrics_corr$treeID)

trees_name1 <- as.character(distances_metrics_corr$treeID)
trees_name2 <- factor(unique(trees_name1))

LAD_metrics1 <- list()
LAD_metrics2 <- list()

for (i in levels(trees_name2)) {
  # Filter data for each tree
  tree1 <- LAD_profiles |> dplyr::filter(treeID == i)
  tree2 <- distances_metrics_corr |> dplyr::filter(treeID == i)

  # Get LAD metrics for each tree
  LAD_metrics <- get_layers_lad(tree1, tree2)
  LAD_metrics1[[i]] <- LAD_metrics$df1
  LAD_metrics2[[i]] <- LAD_metrics$df2
}

LAD_metrics_all1 <- dplyr::bind_rows(LAD_metrics1)
LAD_metrics_all2 <- dplyr::bind_rows(LAD_metrics2)

# List of data frames
LAD_metrics_list <- list(LAD_metrics_all1, LAD_metrics_all2)

# Initialize an empty list to store reordered data frames
reordered_list <- list()

# Specify prefixes (adjust accordingly)
prefixes <- c("treeID", "Hdist", "Hcbh", "effdist", "dptf", "Hdptf", "max", "last")

# Loop over each data frame
for (i in seq_along(LAD_metrics_list)) {

  LAD_metrics_all <- LAD_metrics_list[[i]]

  # Get original column names
  original_column_names <- colnames(LAD_metrics_all)

  # Initialize vector to store new order
  new_order <- c()

  # Loop over prefixes
  for (prefix in prefixes) {
    # Find column names matching the current prefix
    matching_columns <- grep(paste0("^", prefix), original_column_names, value = TRUE)

    # Extract numeric suffixes and order the columns based on these suffixes
    numeric_suffixes <- as.numeric(gsub(paste0("^", prefix), "", matching_columns))
    
       # Order the columns based on numeric suffixes
    matching_columns <- matching_columns[order(numeric_suffixes)]

    # Append to new order
    new_order <- c(new_order, matching_columns)
  }
  # Reorder columns
  LAD_metrics_all <- LAD_metrics_all[, new_order]
  # Store the reordered data frame in the list
  reordered_list[[i]] <- LAD_metrics_all
 }
  # Write the reordered data frame to a file
 fuels_lad_all_path <- file.path(system.file("extdata", package = "LadderFuelsR"), "7_fuels_lad_all.txt")
 write.table(reordered_list[[1]], file = fuels_lad_all_path, sep = "\t", row.names = FALSE)
 fuels_lad_gt25_path <- file.path(system.file("extdata", package = "LadderFuelsR"), "7_fuels_lad_gt25perc.txt")
 write.table(reordered_list[[2]], file = fuels_lad_gt25_path, sep = "\t", row.names = FALSE)
 
 head(reordered_list[[2]])
##                treeID1                treeID Hdist1 Hdist2 Hcbh1 Hcbh2
## Percentage...1       1   1_Eglin_zone1_CROWN    5.5     NA   6.5    NA
## Percentage...2      10  10_Eglin_zone1_CROWN    4.5     NA   5.5    NA
## Percentage...3     100 100_Eglin_zone1_CROWN    4.5     NA   5.5    NA
## Percentage...4     101 101_Eglin_zone1_CROWN    3.5     NA   4.5    NA
## Percentage...5     102 102_Eglin_zone1_CROWN    4.5     NA   5.5    NA
## Percentage...6     103 103_Eglin_zone1_CROWN    4.5     NA   5.5    NA
##                Hcbh1_Hdptf1 Hcbh2_Hdptf2 effdist1 dptf1 dptf2 Hdptf1 Hdptf2
## Percentage...1     94.49448           NA        5     8    NA   14.5     NA
## Percentage...2     99.37928           NA        4     8    NA   13.5     NA
## Percentage...3     94.88224           NA        4     6    NA   11.5     NA
## Percentage...4     99.66928           NA        3     7    NA   11.5     NA
## Percentage...5     98.91892           NA        4     6    NA   11.5     NA
## Percentage...6     95.17946           NA        4     6    NA   11.5     NA
##                max_Hcbh max_dptf max_Hdptf max_height maxlad_Hcbh maxlad_Hdist
## Percentage...1      6.5        8      14.5       15.5         6.5          5.5
## Percentage...2      5.5        8      13.5       14.5         5.5          4.5
## Percentage...3      5.5        6      11.5       12.5         5.5          4.5
## Percentage...4      4.5        7      11.5       12.5         4.5          3.5
## Percentage...5      5.5        6      11.5       12.5         5.5          4.5
## Percentage...6      5.5        6      11.5       12.5         5.5          4.5
##                maxlad_Hdptf maxlad_dptf maxlad_effdist last_Hcbh last_dptf
## Percentage...1         14.5           8              5       6.5         8
## Percentage...2         13.5           8              4       5.5         8
## Percentage...3         11.5           6              4       5.5         6
## Percentage...4         11.5           7              3       4.5         7
## Percentage...5         11.5           6              4       5.5         6
## Percentage...6         11.5           6              4       5.5         6
##                last_Hdptf
## Percentage...1       14.5
## Percentage...2       13.5
## Percentage...3       11.5
## Percentage...4       11.5
## Percentage...5       11.5
## Percentage...6       11.5

#18. PLOT FUEL LAYERS WITH LAD > 25 % AND CBH BASED ON MAXIMUM LAD PERCENTAGE

library(ggplot2)

# LAD profiles derived from normalized ALS data after applying [lad.profile()] function
data_path <- file.path(system.file("extdata", package = "LadderFuelsR"), "LAD_profiles.txt")
LAD_profiles <- read.table(data_path, sep = "\t", header = TRUE)
LAD_profiles$treeID <- factor(LAD_profiles$treeID)

# Tree metrics derived from get_layers_lad() function
LAD_gt25p_path <- file.path(system.file("extdata", package = "LadderFuelsR"), "7_fuels_lad_gt25perc.txt")
fuels_LAD_metrics <- read.table(LAD_gt25p_path, sep = "\t", header = TRUE)
fuels_LAD_metrics$treeID <- factor(fuels_LAD_metrics$treeID)

trees_name1 <- as.character(fuels_LAD_metrics$treeID)
trees_name2 <- factor(unique(trees_name1))

# Generate plots for fuels LAD metrics
plots_trees_LAD <- get_plots_cbh_LAD(LAD_profiles, fuels_LAD_metrics)

# Save plots for each tree
for (name in names(plots_trees_LAD)) {
  plots <- plots_trees_LAD[[name]]

  if (!is.null(plots)) {
    #print(paste("Saving plot for tree:", name))
    ggsave(file.path(system.file("extdata", package = "LadderFuelsR"),  paste0( name, "_LAD_25perc", ".tiff")), plot = plots,
           width = 6, height = 6, units = "in")
  }}

par(mfrow = c(3, 1))
# Plot the first panel
plot(plots_trees_LAD[[1]],width = 6, height = 6)

# Plot the second panel
plot(plots_trees_LAD[[2]],width = 6, height = 6)

# Plot the third panel
plot(plots_trees_LAD[[3]],width = 6, height = 6)

#19. CBH BASED ON THE BREAKING POINT METHOD AND LAD PERCENTAGE

library(dplyr)
library(magrittr)

# LAD profiles derived from normalized ALS data after applying [lad.profile()] function
data_path <- file.path(system.file("extdata", package = "LadderFuelsR"), "LAD_profiles.txt")
LAD_profiles <- read.table(data_path, sep = "\t", header = TRUE)
LAD_profiles$treeID <- factor(LAD_profiles$treeID)

# Tree metrics derived from get_effective_gap() function
distcor_path <- file.path(system.file("extdata", package = "LadderFuelsR"), "6_tree_distances_metrics_corr.txt")
distances_metrics_corr <- read.table(distcor_path, sep = "\t", header = TRUE)
distances_metrics_corr$treeID <- factor(distances_metrics_corr$treeID)

trees_name1 <- as.character(distances_metrics_corr$treeID)
trees_name2 <- factor(unique(trees_name1))

cum_LAD_metrics_list <- list()

for (i in levels(trees_name2)) {
  # Filter data for each tree
  tree1 <- LAD_profiles |> dplyr::filter(treeID == i)
  tree2 <- distances_metrics_corr |> dplyr::filter(treeID == i)

  # Get cumulative LAD metrics for each tree
  cum_LAD_metrics <- get_cum_break(tree1, tree2)
  cum_LAD_metrics_list[[i]] <- cum_LAD_metrics
}
## breakpoint estimate(s): 3.118038
# Combine the individual data frames
cum_LAD_metrics_all <- dplyr::bind_rows(cum_LAD_metrics_list) 

# =======================================================================#
# REORDER COLUMNS
# =======================================================================#

# Get original column names
original_column_names <- colnames(cum_LAD_metrics_all)

# Specify prefixes (adjust accordingly)
prefixes <- c("treeID", "Hcbh", "below", "above", "depth", "Hdepth", "dptf", "Hdptf", "Hdist", "effdist", "max", "last", "cumlad")

# Initialize vector to store new order
new_order <- c()

# Loop over prefixes
for (prefix in prefixes) {
  # Find column names matching the current prefix
  matching_columns <- grep(paste0("^", prefix), original_column_names, value = TRUE)

  # Extract numeric suffixes and order the columns based on these suffixes
  numeric_suffixes <- as.numeric(gsub(paste0("^", prefix), "", matching_columns))
  matching_columns <- matching_columns[order(numeric_suffixes)]

  # Append to new order
  new_order <- c(new_order, matching_columns)
}

# Reorder columns
cum_LAD_metrics_order <- cum_LAD_metrics_all[, new_order]
cum_LAD_metrics_path <- file.path(system.file("extdata", package = "LadderFuelsR"), "8_cbh_breaking_point_lad.txt")
write.table(cum_LAD_metrics_order, file = cum_LAD_metrics_path, sep = "\t", row.names = FALSE)

head(cum_LAD_metrics_order)
##   treeID1                treeID    Hcbh1  Hcbh_bp Hcbh1_Hdptf1 below_hcbhbp
## 1       1   1_Eglin_zone1_CROWN 8.487859 8.487859     89.45101    33.925914
## 2      10  10_Eglin_zone1_CROWN 6.809724 6.809724     80.49151    26.699603
## 3     100 100_Eglin_zone1_CROWN 6.843101 6.843101     96.20392     9.172215
## 4     101 101_Eglin_zone1_CROWN       NA 7.497583           NA    58.726259
## 5     102 102_Eglin_zone1_CROWN 6.367170 6.367170     81.03784    25.827027
## 6     103 103_Eglin_zone1_CROWN 6.168868 6.168868     92.15937    13.009641
##   above_hcbhbp dptf1 Hdptf1 Hdist1 effdist1 maxlad_Hcbh maxlad_Hdptf
## 1     89.45101     7   15.5    7.5        7         8.5         15.5
## 2     80.49151     8   14.5    5.8        6         6.8         14.5
## 3     96.20392     6   12.5    5.8        6         6.8         12.5
## 4     65.31229    NA     NA     NA       NA          NA           NA
## 5     81.03784     6   12.5    5.4        5         6.4         12.5
## 6     92.15937     6   12.5    5.2        5         6.2         12.5
##   maxlad_Hdist maxlad_effdist max_dptf max_Hcbh max_Hdptf max_Hdist max_height
## 1          7.5              7        7      8.5      15.5       7.5       15.5
## 2          5.8              6        8      6.8      14.5       5.8       14.5
## 3          5.8              6        6      6.8      12.5       5.8       12.5
## 4           NA             NA       NA       NA        NA        NA       12.5
## 5          5.4              5        6      6.4      12.5       5.4       12.5
## 6          5.2              5        6      6.2      12.5       5.2       12.5
##   last_effdist last_Hcbh last_Hdptf last_Hdist cumlad
## 1            7       8.5       15.5        7.5 0.4171
## 2            6       6.8       14.5        5.8 0.4013
## 3            6       6.8       12.5        5.8 0.0923
## 4           NA        NA         NA         NA 0.4404
## 5            5       6.4       12.5        5.4 0.2389
## 6            5       6.2       12.5        5.2 0.1120

#20. PLOT CBH BASED ON THE BREAKING POINT METHOD AND LAD PERCENTAGE

library(ggplot2)
# LAD profiles derived from normalized ALS data after applying [lad.profile()] function
data_path <- file.path(system.file("extdata", package = "LadderFuelsR"), "LAD_profiles.txt")
LAD_profiles <- read.table(data_path, sep = "\t", header = TRUE)
LAD_profiles$treeID <- factor(LAD_profiles$treeID)

# Tree metrics derived from get_cum_break() function
cum_LAD_metrics_path <- file.path(system.file("extdata", package = "LadderFuelsR"), "8_cbh_breaking_point_lad.txt")
cum_LAD_metrics <- read.table(cum_LAD_metrics_path, sep = "\t", header = TRUE)
cum_LAD_metrics$treeID <- factor(cum_LAD_metrics$treeID)

# Generate cumulative LAD plots
plots_trees_cumlad <- get_plots_cumm(LAD_profiles, cum_LAD_metrics)

# Save plots for each tree
for (name in names(plots_trees_cumlad)) {
  plots <- plots_trees_cumlad[[name]]

  if (!is.null(plots)) {
    #print(paste("Saving plot for tree:", name))
    ggsave(file.path(system.file("extdata", package = "LadderFuelsR"),  paste0( name, "_cumm_Hcbh", ".tiff")), plot = plots)
  }}

par(mfrow = c(3, 1))
# Plot the first panel
plot(plots_trees_cumlad[[1]],width = 6, height = 6)

# Plot the second panel
plot(plots_trees_cumlad[[2]],width = 6, height = 6)

# Plot the third panel
plot(plots_trees_cumlad[[3]],width = 6, height = 6)

#21. JOINING FUEL LADDER PROPERITES WITH CROWN POLYGONS

# Tree metrics derived from get_layers_lad() function
LAD_gt25p_path <- file.path(system.file("extdata", package = "LadderFuelsR"), "7_fuels_lad_gt25perc.txt")
fuels_LAD_metrics <- read.table(LAD_gt25p_path, sep = "\t", header = TRUE)
fuels_LAD_metrics$treeID <- factor(fuels_LAD_metrics$treeID)

# No overlapping crown polygins (output from step 6)
trees_path <- file.path(system.file("extdata", package = "LadderFuelsR"), "crown2m_silva_convex_no_overlap.shp")
tree_file <- st_read(trees_path)
## Reading layer `crown2m_silva_convex_no_overlap' from data source 
##   `D:\OLGA\R_library\LadderFuelsR\extdata\crown2m_silva_convex_no_overlap.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 877 features and 16 fields
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: 525912.5 ymin: 3378523 xmax: 526028 ymax: 3378614
## Projected CRS: NAD83 / UTM zone 16N
tree_file$treeID1 <- factor(tree_file$treeID)

crowns_properties<-merge (tree_file,fuels_LAD_metrics, by="treeID1")


# Define a red-to-green color palette
red_to_green_palette <- colorRampPalette(c("red", "green"))

# Plotting with the red-to-green color palette
ggplot() +
  geom_sf(data = crowns_properties, aes(fill = Hcbh1)) +
  scale_fill_gradientn(colors = red_to_green_palette(100), na.value = "grey50") +
  theme_minimal() +
  labs(title = "Tree Crowns", fill = "Hcbh1")

crowns_properties$maxlad_Hcbh_factor <- cut(crowns_properties$maxlad_Hcbh, breaks = 5)

# Plotting with a discrete legend
ggplot() +
  geom_sf(data = crowns_properties, aes(fill = maxlad_Hcbh_factor)) +
  scale_fill_manual(values = red_to_green_palette(5)) +
  theme_minimal() +
  labs(title = "Tree Crowns", fill = "maxlad_Hcbh")